package scu.maqiang.bes;

import java.util.Arrays;
import java.util.HashSet;

import scu.maqiang.mesh.Mesh2T3;
import scu.maqiang.numeric.*;

public class BES2DT3L20 extends BES {
	public BES2DT3L20(Mesh2T3 mesh2D) {
		super(mesh2D, 1);
		xM = new double[nB];
		yM = new double[nB];
		for (int i = 0; i < nB; i++) {
			int[] ele = mesh.getBoundary(i);
			mesh.getCoordPerBoundary(ele, coordBoundary);
			xM[i] = 0.5 * (coordBoundary[0][0] + coordBoundary[1][0]);
			yM[i] = 0.5 * (coordBoundary[0][1] + coordBoundary[1][1]);
		}
		mesh.getBoundaryStream().forEach(ele -> {
			for (int i = 0, n = ele.length; i < n; i++) {
				mesh.setNodeLabel(i, 10);
			}
		});
		knownValue = new double[ndof];
	}

	public int GetNdof() {
		return ndof;
	}

	public void setDirichletBoundaryValue(double value, int... dirichletBoundary) {
		for (int i = 0, nd = dirichletBoundary.length; i < nd; i++) {
			dirichletBd.add(dirichletBoundary[i]);
		}

		for (int i = 0; i < nB; i++) {
			int label = mesh.getBoundaryLabel(i);
			for (int j = 0, nd = dirichletBoundary.length; j < nd; j++) {
				if (label == dirichletBoundary[j]) {
					knownValue[i] = value;
					break;
				}
			}
		}
	}

	public void setNeumannBoundaryValue(double value, int... neumannBoundary) {
		for (int i = 0, nd = neumannBoundary.length; i < nd; i++) {
			neumannBd.add(neumannBoundary[i]);
		}

		for (int i = 0; i < nB; i++) {
			int label = mesh.getBoundaryLabel(i);
			for (int j = 0, nd = neumannBoundary.length; j < nd; j++) {
				if (label == neumannBoundary[j]) {
					knownValue[i] = value;
					break;
				}
			}
		}
	}

	public double[] getKnownValue() {
		return knownValue;
	}

	public void assembleSource(double[] constSource, double[] source) {
		double[][] point = Constants.TP5;
		double[] weight = Constants.TW5;
		Arrays.fill(source, 0.0);
		//Mesh2DT3Generator mg = new Mesh2DT3Generator();
		//Mesh2DT3 mesh = mg.square2D(4, 4, T3Type.Left);
		//mesh.displayNodes();
		//mesh.displayElements();
		//System.exit(0);
		for(int i = 0; i < nB; i++) {
			double xp = xM[i], yp = yM[i];
			for(int j = 0, n = mesh.getNt(); j < n; j++) {
				mesh.getCoordPerElement(mesh.getElement(j), coordRegion);
				double area = CG.area2DT3(coordRegion);
				double[] x = MVO.col(coordRegion, 0);
				double[] y = MVO.col(coordRegion, 1);
				for(int k = 0, nn = weight.length; k < nn; k++) {
					double xK = MVO.dot_product(x, point[k]);
					double yK = MVO.dot_product(y, point[k]);
					double r = Math.sqrt((xK - xp) * (xK - xp) + (yK - yp) * (yK - yp));
					//System.out.println(xp + "\t" + yp + "\t" + xK + "\t" + yK + "\t" + area + "\t" + r);
					source[i] += -area * Math.log(r) * constSource[0] * weight[k] / (2 * Math.PI);
				}
				//System.out.println();
			}
		}
		
	}
	public void assembleLaplaceMatrix(double[] constCoef, double[][] H, double[][] G) {
		double[][] GaussQuadrature = Constants.GaussQuadratureConstants[3];
		int order = GaussQuadrature[0].length;
		double[] GP = GaussQuadrature[0];
		double[] GW = GaussQuadrature[1];
		double H_ij = 0.0, G_ij = 0.0;
		double DistIJ = 0.0;
		for (int i = 0; i < nB; i++) {
			double xm = xM[i], ym = yM[i];
			for (int j = 0; j < nB; j++) {
				mesh.getCoordPerBoundary(mesh.getBoundary(j), coordBoundary);
				double length = CG.length2DL2(coordBoundary);
				double[] x = { coordBoundary[0][0], coordBoundary[1][0] };
				double[] y = { coordBoundary[0][1], coordBoundary[1][1] };
				if (Math.abs(x[0] - x[1]) < 1.0e-10) {
					DistIJ = Math.abs(xm - x[0]);
				} else {
					DistIJ = Math.abs((x[0] - xm) * (y[1] - y[0]) - (y[0] - ym) * (x[1] - x[0])) / length;
				}
				DistIJ = CG.area2DT3(new double[] { xm, ym }, coordBoundary[0], coordBoundary[1]) > 0 ? DistIJ : -DistIJ;
				if (i == j) {
					H_ij = 0.5;
					G_ij = length / (2 * Math.PI) * (1 - Math.log(0.5 * length));
				} else {
					H_ij = 0.0;
					G_ij = 0.0;
					for (int k = 0; k < order; k++) {
						double x0 = 0.5 * (1 - GP[k]) * x[0] + 0.5 * (1 + GP[k]) * x[1];
						double y0 = 0.5 * (1 - GP[k]) * y[0] + 0.5 * (1 + GP[k]) * y[1];
						double rij = Math.sqrt((xm - x0) * (xm - x0) + (ym - y0) * (ym - y0));
						// System.out.println(Arrays.toString(GP));
						//System.out.println(xm + "\t" + ym + "\t" + x0 + "\t" + y0 + "\t" + rij);
						// System.exit(0);
						H_ij += GW[k] / (rij * rij);
						G_ij += GW[k] * Math.log(rij);
					}
					H_ij *= -0.25 * length / Math.PI * DistIJ;
					G_ij *= -0.25 * length / Math.PI;
				}
				//System.out.println(i + "\t" + j + "\t" + H_ij + "\t" + G_ij);
				if (dirichletBd.contains(mesh.getBoundaryLabel(j))) {
					H[i][j] = -G_ij;
					G[i][j] = -H_ij;
				} else {
					H[i][j] = H_ij;
					G[i][j] = G_ij;
				}
			}
		}
	}

	public double[][] getSolution(double[] sol) {
		double[][] result = new double[2][ndof];
		for (int i = 0; i < nB; i++) {
			if (dirichletBd.contains(mesh.getBoundaryLabel(i))) {
				result[0][i] = knownValue[i];
				result[1][i] = sol[i];
			} else {
				result[0][i] = sol[i];
				result[1][i] = knownValue[i];
			}
		}
		return result;
	}

	public void displayMiddlePoint() {
		for (int i = 0; i < nB; i++) {
			System.out.println(xM[i] + "\t" + yM[i]);
		}
	}

	public double computeValue(double[] sol, double[] dsol, double[] source, double[] xy) {
		double[][] GaussQuadrature = Constants.GaussQuadratureConstants[3];
		int order = GaussQuadrature[0].length;
		double[] GP = GaussQuadrature[0];
		double[] GW = GaussQuadrature[1];
		double H_ij = 0.0, G_ij = 0.0;
		double DistIJ = 0.0;
		double result = 0.0;
		double xm = xy[0], ym = xy[1];
		double[][] point = Constants.TP5;
		double[] weight = Constants.TW5;
		for (int j = 0; j < nB; j++) {
			mesh.getCoordPerBoundary(mesh.getBoundary(j), coordBoundary);
			double length = CG.length2DL2(coordBoundary);
			double[] x = { coordBoundary[0][0], coordBoundary[1][0] };
			double[] y = { coordBoundary[0][1], coordBoundary[1][1] };
			if (Math.abs(x[0] - x[1]) < 1.0e-10) {
				DistIJ = Math.abs(xm - x[0]);
			} else {
				DistIJ = Math.abs((x[0] - xm) * (y[1] - y[0]) - (y[0] - ym) * (x[1] - x[0])) / length;
			}
			DistIJ = CG.area2DT3(xy, coordBoundary[0], coordBoundary[1]) > 0 ? DistIJ : -DistIJ;
			H_ij = 0.0;
			G_ij = 0.0;
			for (int kk = 0; kk < order; kk++) {
				double x0 = 0.5 * (1 - GP[kk]) * x[0] + 0.5 * (1 + GP[kk]) * x[1];
				double y0 = 0.5 * (1 - GP[kk]) * y[0] + 0.5 * (1 + GP[kk]) * y[1];
				double rij = Math.sqrt((xm - x0) * (xm - x0) + (ym - y0) * (ym - y0));
				H_ij += GW[kk] / (rij * rij);
				G_ij += GW[kk] * Math.log(rij);
			}
			H_ij *= -0.25 * length / Math.PI * DistIJ;
			G_ij *= -0.25 * length / Math.PI;
			result += dsol[j] * G_ij - sol[j] * H_ij;
		}
		if (source != null) {
			double BI = 0.0;
			//Mesh2DT3Generator mg = new Mesh2DT3Generator();
			//Mesh2DT3 mesh = mg.square2D(4, 4, T3Type.Left);
			for(int jj = 0, n = mesh.getNt(); jj < n; jj++) {
				mesh.getCoordPerElement(mesh.getElement(jj), coordRegion);
				double area = CG.area2DT3(coordRegion);
				double[] x = MVO.col(coordRegion, 0);
				double[] y = MVO.col(coordRegion, 1);
				for(int k = 0, nn = weight.length; k < nn; k++) {
					double xK = MVO.dot_product(x, point[k]);
					double yK = MVO.dot_product(y, point[k]);
					double r = Math.sqrt((xK - xm) * (xK - xm) + (yK - ym) * (yK - ym));
						//System.out.println(xp + "\t" + yp + "\t" + xK + "\t" + yK + "\t" + area + "\t" + r);
					BI += -area * Math.log(r) * source[0] * weight[k] / (2 * Math.PI);
				}
					//System.out.println();
			}
			result += BI;
			
		}
		return result;
	}

	public double[][] computeDomainValues(double[] sol, double[] dsol, double[] source) {
		double[][] result = new double[2][mesh.getNv()];
		double[] angles = mesh.computeAngles();
		//System.out.println(MVO.toString(angles));
		for(int i = 0; i < mesh.getNv(); i++) {
			double[] node = mesh.getNode(i);
			result[0][i] = computeValue(sol, dsol, source, node);
			result[0][i] *= 2 * Math.PI / angles[i];
		}
//		for(int i = 0; i < nB; i++) {
//			int ele[] = mesh.getBoundary(i);
//			if(dirichletBd.contains(mesh.getBoundaryLabel(i))) {
//				result[0][ele[0]] = sol[i];
//				result[0][ele[1]] = sol[i];
//			}
//		}
		return result;
	}

	double[] xM = null;
	double[] yM = null;
	HashSet<Integer> dirichletBd = new HashSet<>();
	HashSet<Integer> neumannBd = new HashSet<>();
	double[] knownValue = null;
}
